
#include "Geometry.hpp"
#include <fstream>
using namespace std;
Geometry::Geometry(const std::string &fname)
{
    Massert error;
    fnameGrid = fname;
    ifstream fgrid{fnameGrid};
    error.Assert(fgrid.is_open(), "open grid file error");

    fgrid >> numNode >> numTria;
    coord = new Point[numNode]{};
    elem = new Tria[numTria]{};
    area = new Area[numTria]{};

    for (long i = 0; i < numNode; ++i)
    {
        fgrid >> coord[i].x >> coord[i].y >> coord[i].z;
    }
    for (long i = 0; i < numTria; ++i)
    {
        fgrid >> elem[i].node[0] >> elem[i].node[1] >> elem[i].node[2];
        // 修正为从 0 开始
        elem[i].node[0] -= 1;
        elem[i].node[1] -= 1;
        elem[i].node[2] -= 1;
    }
    fgrid.close();
}

Geometry::~Geometry()
{
    delete[] coord;
    delete[] elem;
    delete[] area;
}

void Geometry::CalcuArea(const Point check)
{
    Point p1, p2, p0, center, vecCheck;
    double tmp;
    for(long i = 0; i < numTria; ++i)
    {
        p0 = coord[elem[i].node[0]];
        p1 = coord[elem[i].node[1]];
        p2 = coord[elem[i].node[2]];

        p1.x -= p0.x;
        p1.y -= p0.y;
        p1.z -= p0.z;
        p2.x -= p0.x;
        p2.y -= p0.y;
        p2.z -= p0.z;

        center.x = (p0.x + p1.x + p2.x)/3.0;
        center.y = (p0.y + p1.y + p2.y)/3.0;
        center.z = (p0.z + p1.z + p2.z)/3.0;

        vecCheck.x = center.x - check.x;
        vecCheck.y = center.y - check.y;
        vecCheck.z = center.z - check.z;

        area[i].x = 0.5*(p1.y*p2.z - p1.z*p2.y);
        area[i].y = 0.5*(p1.z*p2.x - p1.x*p2.z);
        area[i].z = 0.5*(p1.x*p2.y - p1.y*p2.x);

        tmp = area[i].x*vecCheck.x + area[i].y*vecCheck.y + area[i].z*vecCheck.z;
        if(tmp < 0.0)
        {
            area[i].x = -area[i].x;
            area[i].y = -area[i].y;
            area[i].z = -area[i].z;
        }
    }
}